function [D X Di Xi] = ResultsEP2(output, subjectNrs, electrode, period)
X = []; D = [];
Xi = []; Di = [];
% skipnext = true;
smoothCount = 0;
%% Normalizing per subject (each period for a certain electrode)
% loop over subjects which checks if the current period exists in the
% current subject and if so runs ResultsSEP on it, which creates a
% 'data' matrix from the output/results matrix structure, normalizes this new matrix (data), and smoothes the normalized matrix (so it smoothes per subject).
for i=subjectNrs
    disp(['Subject: ' output.subjects(i).Name]);
    %     if skipnext; continue; end;
    
    if length(output.results(i).period)<period
        Xep = []; Dep = [];
        disp('     - Period unavailable...');
    else
        [Dep Xep Diep Xiep] = ResultsSEP(output, i, electrode, period, true); % Data per subject is already normalized
        % Diep/Xiep are interpolated variants
        smoothCount = smoothCount + 1;
    end
    D = [D Dep];
    X = [X Xep];
    
    Di(:,:,i) = Diep; % Build the 3d matrix containing all subjects: 8 (lat+val+lat+val...) x 1000 datapoints x number of subjects
    Xi = Xiep;
end
[sy si] = sort(X);
D = D(:,si);
X = X(si);



%% smoothing over all subjects
% maartje update 31/05/10. Remove smoothing over all subjects (leaving only
% the smoothing per subject in ResultsSEP)

% if ~isempty(D)
%     window = 30*smoothCount;
%     if window > (size(X,2)-1)
%         window = (size(X,2)-1);
%     end
%     disp('    Smoothing');
%     for d=1:8
% 
%         % Version with fix for peaks at edges ((r) Roel)
%         temp = [repmat(mean(D(d,1:window)),1,window) D(d,:) repmat(mean(D(d,end-window:end)),1,window)];
%         temp = fastsmooth(temp(:), window, 3, 1);
%         D(d,:) = temp(window+1:end-window);
% 
%         % Without fix
%         %     D(d,:) = fastsmooth(D(d,:), window, 3, 1);
%     end
% end
end